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PREFACE 


The  study  reported  herein  was  sponsored  by  the  Office,  Chief 
of  Engineers,  U.  S.  Army,  as  part  of  the  Civil  Works  General 
Investigations,  Environmental  and  Water  Quality  Operational 
Studies  (EWQOS)  Program.  Work  Unit  No.  31593  (Task  IA.4)  enti¬ 
tled  "Improve  and  Verify  Multidimensional  Hydrodynamic  Mathemat¬ 
ical  Models  for  Reservoirs"  supported  the  subject  study. 

The  study  was  conducted  during  the  period  July  1978  to 
October  1980  by  Dr.  John  E.  Edinger  and  Mr.  Edward  M.  Buchak  of 
J.  E.  Edinger  Associates,  Inc.  Mr.  Mark  S.  Dortch  and  Dr.  Billy 
H.  Johnson  of  the  Hydraulics  Laboratory,  U.  S.  Army  Engineer 
Waterways  Experiment  Station  (WES)  monitored  the  effort.  This 
report  was  written  by  Dr.  Edinger  and  Mr.  Buchak.  Program  man¬ 
ager  of  EWQOS  was  Dr.  Jerome  L.  Mahloch,  WES  Environmental 
Laboratory . 

Commanders  and  Directors  of  WES  during  this  study  and  the 
preparation  of  this  report  were  COL  John  L.  Cannon,  CE, 

COL  Nelson  P.  Conover,  CE,  and  COL  Tilford  C.  Creel,  CE.  Tech¬ 
nical  Director  was  Mr.  Fred  R.  Brown. 


This  report  should  be  cited  as  follows: 

Edinger,  J.  E.,  and  E.  M.  Buchak.  1983.  "Developments 
in  LARM2:  A  Longitudinal-Vertical,  Time-Varying  Hydro- 
dynamic  Reservoir  Model,"  Technical  Report  E-83-1,  pre¬ 
pared  by  J.  E.  Edinger  Associates,  Inc.,  for  the  U.  S. 
Army  Engineer  Waterways  Experiment  Station,  Vicksburg, 
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DEVELOPMENTS  IN  LARM2 :  A  LONGITUDINAL- VERTICAL . 

TIME-VARYING  HYDRODYNAMIC  RESERVOIR  MODEL 

1.  INTRODUCTION 

LARM  developments  for  1979-1980  Include  (1)  modifications 
and  refinements  to  produce  LARM2  as  presented  in  the  User  Guide 
(Buchak  and  Edinger,  1982);  (2)  numerous  investigations  with 
LARM2 , including  sensitivity  analyses  and  testing;  (3)  develop¬ 
ment  of  a  water  quality  transport  module  (WQTM)  for  use  in  mul¬ 
tiple  and  interacting  water  quality  constituent  transport;  (4) 
investigation  of  special  topics  related  to  LARM-type  simulations, 
including  eddy  coefficient  evaluation,  hydrodynamic  volume  con¬ 
veyance,  and  branching;  and  (5)  recommendations  for  further 
extensions. 

The  USER  GUIDE  FOR  LARM2;  A  LONGITUDINAL- VERTICAL,  TIME- 
VARYING  HYDRODYNAMIC  RESERVOIR  MODEL  includes  (1)  the  latest  LARM2 
modifications  for  upstream  cell  addition  and  subtractions  and  (2) 
surface  volume  addition/subtraction  corrections  due  to  depth- 
variable  width.  The  code  and  user  guide  also  include  new  rou¬ 
tines  for  tributary  inflow  and  withdrawal  computations.  The 
present  report  covers  the  remaining  subject  areas  of  LARM2  in¬ 
vestigations,  the  WQTM  with  sediment  and  constituent  transport 
as  examples,  and  the  special  topics  related  to  future  LARM 
development . 
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2.  INVESTIGATIONS  WITH  LARM2 

Numerous  case  studies  have  been  made  with  LARM2,  and  each 
has  provided  new  insight  to  its  applicability.  A  number  of  case 
studies  by  J.E.  Edinger  Associates,  Inc.  (JEEAI)  and  others  are 
on  Corps-related  projects.  The  recent  JEEAI  reservoir  studies 
are  reported  in  the  User  Guide  (Buchak  and  Edinger,  1982).  The 
Corps  related  studies  are  by  Gordon  (1979),  Johnson  (1981),  and 
Edinger  and  Buchak  (1980).  These  studies  have  shown  the  versa¬ 
tility  of  the  LARM2  computational  basis  and  have  indicated  para¬ 
meter  sensitivity  under  a  wide  range  of  conditions. 

The  studies  by  Gordon  (1979)  of  the  Center  Hill  Reservoir 
include  a  sensitivity  study  as  well  as  a  seasonal  temperature 
verification.  The  investigations  by  Johnson  (1981)  included 
simulation  of  an  underflow  density  current  in  a  Waterways 
Experiment  Station  (WES)  flume  and  determined  the  parameter 
effects  of  scaling  to  small  grid  sizes  of  transitional  laminar 
flows.  The  studies  by  Edinger  and  Buchak  (1980)  include  the 
conversion  of  LARM2  reservoir  boundary  conditions  to  estuarine 
boundaries  and  the  comparison  of  computed  and  observed  veloc¬ 
ities  for  tidal  flows.  Each  of  these  studies  has  shown  the 
sensitivity  of  LARM  computations  to  different  parameters  under 
different  conditions. 

Center  Hill  Studies 

Center  Hill  Lake  near  Cookeville,  Tennessee,  was  studied  by 
Gordon  (1979)  using  LARM2.  One  purpose  of  the  study  was  to  verify 


the  ability  of  LARM2  to  predict  observed  temperature  distribu¬ 
tions  using  only  geometry,  inflow,  and  meteorological  data. 
Another  purpose  of  the  study  was  to  determine  the  sensitivity 
of  LARM2  velocity  profiles  to  different  types  of  outfall 
conditions . 

The  temperature  verifications  of  LARM2,  as  in  other  cases 
were  quite  successful.  A  comparison  of  observed  and  predicted 
temperatures  for  the  Center  Hill  Lake  simulations  are  shown  in 
Tables  2.1  to  2.4.  Some  sensitivity  was  shown  to  input  data 
anomalies,  as  would  be  expected.  There  were  no  measurements  of 
velocity  profiles  in  the  reservoir,  and  it  can  only  be  assumed 
that  overall  circulation  features  were  correct  since  the  heat 
budget  was  maintained. 

Center  Hill  Lake  has  a  low-level  outlet  about  30  m  below 
the  water  surface  that  normally  withdraws  hypolimnetic  water. 
Simulations  were  made  with  an  outlet  at  10  m  below  the  water 
surface  to  withdraw  metalimnetic  water,  and  for  a  hypothetical 
submerged  weir  upstream  of  the  dam. 

For  the  case  of  the  high-level  outlet  Gordon  (1979)  con¬ 
cludes  from  LARM2  output: 

(1)  The  hypolimnion  stays  much  cooler  throughout 
the  simulation  period,  as  overlying  water  is 
discharged.  (Both  normal  and  raised-outlet 
simulations  were  started  with  identical 
stratification  conditions.) 

(2)  Discharge  temperatures  are  considerably 
warmer  with  raised  outlet. 


(3)  Inflow  patterns  continued  to  be  the  same 
as  for  the  low-level  outlet,  with  surface 
flows  in  the  spring  and  shallow  inter¬ 
flows  during  the  summer. 

Although  certain  of  these  conclusions  might  have  been  inferred 
without  LARM,  their  demonstration  is  quite  striking. 

The  submerged  weir  simulation  was  made  by  "zeroing  out" 
horizontal  velocities  and  dispersion  terms  along  the  interface 
between  cells  extending  from  the  bottom  to  within  10  m  of  the 
surface.  The  algorithm  for  this  change  was  quite  straight¬ 
forward  and  minimal.  The  longitudinal-vertical  vector  plots 
indicated  velocities  that  looked  like  flows  over  a  submerged 
weir. 

Gordon  (1979)  expected  to  find  temperpture  discontinuities 
on  either  side  of  the  weir  as  the  downstream  volume  filled  with 
surface  water  over  the  weir.  It  was  found,  however,  that  since 
inflows  were  either  at  the  surface  or  shallow  underflows,  the 
submerged  weir  had  little  effect  on  downstream  temperatures. 


Potomac  Estuary  Studies 

The  LARM2  hydrodynamic  boundary  conditions  were  modified 
to  compute  as  an  inflow,  downestuary  head  problem  with  the  down- 
estuary  head  being  a  specified  tide.  At  the  tidal  boundary,  the 
vertical  profile  of  velocity  is  computed  using  specified  water 
surface  elevations  and  a  specified  salinity  profile.  Salinity 
transport  was  computed  in  place  of  heat  transport.  The 


estuarine  version  of  LARM2  is  called  LAEM  (Laterally  Averaged 
Estuary  Model)  and  has  been  tested  for  the  Potomac  River  estu¬ 
ary  (Edinger  and  Buchak,  1980). 

The  estuary  simulations  are  important  because  estuary 
velocities  reach  easily  detectable  levels  (on  the  order  of 
30  cm/s),  and  computed  velocities  can  be  compared  to  measured 
velocities.  Also,  water  surface  elevations  change  rapidly  with 
time,  and  significant  longitudinal  water  surface  profiles 
develop.  Reservoirs  have  low  velocities  and  water  surfaces 
that  change  mostly  with  storage.  Although  LARM2  has  been  tested 
on  numerous  occasions  against  measured  temperature  fields  in 
reservoirs,  these  simulations  have  not  provided  a  rigorous  test 
of  dynamic  and  velocity  computations.  Other  features  of  strat¬ 
ified  estuaries,  such  as  the  distribution  of  vertical  diffusiv- 
ities  at  dif:  srent  parts  of  the  tidal  cycle  are  well  known  and 
can  be  compared  to  computed  values. 

The  Potomac  River  estuary  extends  167  km  from  Chesapeake 
Bay  to  Great  Falls,  upestuary  of  Washington,  DC.  The  "ocean" 
relative  to  the  estuary  is  Chesapeake  Bay,  and  the  estuary  re¬ 
sponds  to  the  tide  height  variation  and  vertical  salinity  dis¬ 
tribution  in  the  bay.  During  September  1974,  hourly  velocity 
and  salinity  profiles  were  measured  at  stations  that  were 
19  and  35  km  from  the  estuary  mouth  (Stations  P10  and  P19, 
respectively).  Data  at  these  stations  are  used  for  verifying 
the  ability  of  LAEM  to  compute  velocity  fields. 


The  LAEM  simulation  was  structured  with  a  Ax  of  9  km  and 
layer  thicknesses  of  1  m.  A  single-period  tide  of  0.2  m  amplitude 
and  salinity  profile  were  specified  at  the  mouth.  September  1974 
freshwater  inflows  were  specified  at  the  upestuary  end. 

The  narrowing  cross-sectional  geometry  of  the  Potomac  is 
such  that  the  tidal  amplitude  generally  increases  in  the  up¬ 
estuary  direction.  A  rapid  change  in  geometry  about  93  km  from 
the  mouth  reduces  the  tidal  amplitude  to  its  minimum.  The  geom¬ 
etry  of  the  estuary  also  modifies  the  phase  lag  as  the  tidal 
wave  progresses  up  the  estuary. 

The  first  two  properties  of  the  estuary  compared  to  the 
hydrodynamic  computations  are  the  observed  tidal  range  and 
tidal  phase  lag.  (The  observed  tidal  ranges  presented  in 
Figure  2.1  are  the  maximum  values  that  occurred  independently 
at  each  station  over  a  number  of  years  and  are  not  values  con¬ 
current  with  the  predictions.  )  The  comparisons  are  shown  in 
Figures  2.1  and  2.2.  The  computations  are  presented  for  dif¬ 
ferent  values  of  the  bottom  friction  Chezy  coefficient.  The 
latter  is  taken  as  a  constant  over  the  length  of  the  estuary,  al¬ 
though  the  LARM  code  does  provide  for  spatially  varying  friction 
coefficients.  It  is  found  that  a  relatively  high  Chezy  coeffi- 
cent  is  generally  required  to  reproduce  the  observed  longitudinal 
variation  in  tidal  range. 

The  observed  and  computed  phase  lags  are  compared  in 
Figure  2.2.  The  phase  lag  is  less  sensitive  to  friction  than  is 


tidal  amplitude  because  phase  lag  is  mostly  a  function  of  tidal 
wave  speed  and  hence  depth.  The  comparison  of  observed  and  com¬ 
puted  values  is  an  indication  of  how  well  the  LAEM  geometry 
schematization  represents  the  real  waterbody  dynamically. 

Comparisons  of  observed  (+)  and  computed  velocities  at 
four  depths  for  the  station  35  km  from  the  mouth  are  given  in 
Figure  2.3  for  two  tidal  cycles.  The  simulations  were  ini¬ 
tialized  over  a  period  of  ten  tidal  cycles.  The  simulations 
reproduce  the  velocity  amplitudes  and  their  decrease  with  depth. 
There  are  some  phase  shifts  between  observed  and  computed 
velocities  due  to  short  term  wind  effects  not  included  in  the 
model  tide  representation  at  the  mouth  of  the  estuary.  The 
results  show  that  the  LAEM  dynamics  adequately  reproduce 
rapidly  varying  velocity  conditions. 

A  more  stringent  test  for  a  tidal  estuary  is  reproduction 
of  the  tidally  averaged  velocity.  The  tidally  averaged  velocity 
is  often  an  order  of  magnitude  smaller  than  the  velocity  ampli¬ 
tude  and  approaches  velocities  typically  found  in  reservoirs. 

The  tidally  averaged  velocity  comparisons  are  shown  in 
Figure  2.4  and  show  faithful  reproduction.  The  surface  re¬ 
versal  of  observed  velocities  is  due  to  wind  effects  not  in¬ 
cluded  in  the  simulations.  An  important  feature  of  the  velocity 
profile  is  the  zero  velocity  crossing  point.  It  is  reproduced 
by  the  model. 


The  vertical  distribution  of  the  eddy  viscosity  and  dis¬ 
persion  coefficient  varies  throughout  a  tidal  cycle  due  to 
varying  velocity  shear  relative  to  stationary  salinity  gradient. 
The  computed  dispersion  coefficient  profiles  are  shown  for 
each  tidal  hour  in  Figure  2.5.  The  coefficient  profiles  com¬ 
puted  by  LAEM  vary  throughout  the  cycle  as  expected  from  other 
estuary  studies.  Other  methods  of  formulation  are  considered 
in  Chapter  4. 

Sensitivity  Tests  with  WES  GRH  Flume 

LARM2  was  applied  to  the  WES  Generalized  Reservoir  Hydro¬ 
dynamics  (GRH)  flume  in  order  to  compare  computed  velocity  and 
temperature  fields  with  those  observed  during  several  experi¬ 
ments  conducted  in  1979.  The  well- controlled  conditions  at 
this  facility  and  the  ability  to  obtain  directly  measured  ve¬ 
locity  profiles  made  this  a  unique  application.  These  com¬ 
parisons  resulted  in  sensitivity  tests  for  scaling  of  disper¬ 
sion  parameters,  the  introduction  of  the  Richardson  number 
formulation  for  vertical  dispersion  coefficients,  and  the  use 
of  upstream  differencing  for  the  advection  of  momentum  term  in 
the  longitudinal  momentum  balance. 

The  GRH  flume  and  experiment  are  described  in  Johnson  (1981). 
Basically,  the  experiment  consisted  of  an  upstream  release  of 
16.15*C  water  into  the  flume  which  was  filled  with  21.4 *C  water. 
The  flume  was  initially  at  rest,  and  flow  was  initiated  by 
setting  an  inflow  rate  of  0.00063  m  /s.  The  outflow  was  taken 
from  near  the  bottom  at  the  same  rate  as  the  inflow. 
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The  LARM2  computational  grid  chosen  by  WES  for  simulation 
of  the  experiments  was  Ax*=1.524  m  and  h=0.0762  m.  A  character¬ 
istic  velocity  for  the  flume  is  the  outlet  velocity  of  0.0165  m/s. 
Based  on  a  molecular  kinematic  viscosity  of  water  of  1.5xl0-sm2/s, 
the  Reynolds  number  for  a  computational  cell  is  840,  which  is 
well  within  the  laminar  flow  range.  The  LARM2  simulations  were, 
therefore,  based  on  setting  the  longitudinal  dispersion  coeffi¬ 
cients  of  momentum  and  heat  to  their  molecular  diffusion  values 
of  Ax  =  1.5xio  6  m2/s  and  Dx  =  1.4xio~5  m2/s. 

The  numerical  stability  limits  for  the  LARM2  implicit 
scheme  based  on  the  above  dimensions  are: 

(1)  At  <  Ax/U  -  92  s 

(2)  At  <  Ax  /(2Ax)  -  7.  7  x  10s  s 

(3)  At  <  Ax  /<2Dx)  -  8.3X10“  s 

(4)  At  <  h2/ (2D  )  ■  207  s  for  D  at  its  molecular  value 

z  z 

Clearly,  the  Torrence  condition  given  in  (1)  dominates  the  choice 
of  time  step.  Note  that  the  Courant  or  gravity  wave  speed  crite¬ 
rion  limits  an  explicit  computational  scheme  to  At  <0.51  s. 

The  first  simulation  was  for  A  and  D  at  their  molecular 

X  X 

values,  with  determined  from  a  Richardson  number  criterion  in 
regions  of  unstable  vertical  density  gradients.  The  At  was 
initially  taken  as  60  s.  In  the  initial  time  steps,  the  simu¬ 
lation  immediately  produced  a  characteristic  density  underflow 
at  the  head  of  the  channel  showing  a  velocity  and  temperature 
front  moving  along  the  bottom.  As  the  density  front  velocity 
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increased,  however,  the  Torrence  condition  was  exceeded  locally 
and  computations  at  a  At=60  s  were  terminated.  The  At  was  suc¬ 
cessively  lowered  until  the  maximum  computed  velocity  satisfied 
the  Torrence  condition.  This  was  at  a  At=15  s,  with  a  maximum 
local  velocity  in  the  wedge  of  0.06  m/s. 

The  simulations  at  At=15  s  eventually  produced  local  tem¬ 
perature  inversions  (anomalies  of  0.1  to  0.3  °C) .  It  was  obvious 
that  the  At  was  too  large  for  local  convective  mixing  to  be  com¬ 
pleted  relative  to  advection.  The  theoretical  basis  for  this 
effect  is  found  by  combining  relations  (1)  and  (4)  above.  These 
anomalies  were  eliminated  when  At  was  reduced  to  5  s,  its  final 
value. 

The  second  simulation  set  was  based  on  evaluating  the  shear 
stress  terms  from  the  velocity  gradient  and  evaluating  the  verti¬ 
cal  eddy  viscosity  from  a  Richardson  number  (Ri)  criterion.  For 
high  positive  values  of  the  Ri,  the  lower  limit  to  was  taken 
as  the  molecular  value.  For  high  negative  values,  the  upper 

limit  was  taken  as  A  <  hz/2At  to  be  compatible  with  condition 

z 

(4)  above  and  as  used  for  the  Richardson  criterion  on  vertical 
mixing  in  the  heat  balance.  These  simulations  showed  substan¬ 
tially  the  same  results  as  the  first. 

An  examination  of  the  magnitude  of  the  individual  terms  in 
the  x-direction  momentum  equation  showed  that  the  advection  of 
velocity  terras  (t-—  and  |^)  in  the  previous  simulations  of  the 

dX  dZ 

flume  were  of  the  same  order  of  magnitude  as  the  main  driving 
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3p. 


term,  the  horizontal  pressure  gradient  (^) 


These  simulations 
were  the  first  in  which  the  advection  terms  were  of  such  impor¬ 
tance.  The  velocity  fields  computed  in  the  first  two  simulation 
sets  also  showed  the  formation  of  prominent  eddies,  which  were 


not  apparent  in  the  flume  test.  For  these  two  reasons,  it  was 
decided  to  apply  upstream  differencing  throughout  the  grid  for 
both  the  -|^  and  -|^  terms  instead  of  only  for  the  term  in 
the  vicinity  of  the  outlet.  This  change  resulted  in  an  improved 
time-temperature  curve  and  more  rational  velocity  fields. 


Km  from  Head  of  Tide 


Figure  2.1 

Comparison  of  Observed  and 
Predicted  Tidal  Ranges  along 
the  Potomac  Estuary 


Km  from  Head  of  Tide 


Figure  2.2 

Comparison  of  Observed  and 
Predicted  Tidal  Phase  along 
the  Potomac  Estuary 


Figure  2 . 3 

Comparison  of  Observed  (+)  and  Computed 
Velocities  35  km  from  the  Estuary 
Mouth  (Station  P19)  for  the 
Indicated  Depths 


TABLE  2.1  Center  Hill  Lake  Comparison  of  Predicted  and  Observed  Temperatures 
(°C)  for  April  4,  1972  (Gordon,  1979). 


(2)  Type  of  Simulation  -  Isothermal  to  spring  stratification 

(3)  Other  -  Average  mean  error  of  simulation  ■  -0.54°C 


TABLE 


(2)  Type  of  Simulation  -  Summer  stratification 

(3)  Other  -  Average  mean  error  of  simulation  ■  +1.13°C 


TABLE  2.3  Center  Hill  Lake  Comparison  of  Predicted  and  Observed  Temperatures 
( °C)  for  March  29,  1973  (Gordon,  1979). 


TABLE  2.4  Center  Hill  Lake  Comparison  of  Predicted  and  Observed  Temperat 
(°C)  for  August  15,  1973  (Gordon,  1979). 


(1)  Days  into  run  sequence  -  51 

(2)  Type  of  Simulation  -  Samer  stratification 

(3)  Other  -  Average  aean  eirot  of  simulation  ■  +0.73°C 


3.  THE  WATER  QUALITY  TRANSPORT  MODULE 
AND  CONSTITUENT  SOURCE-SINK  ROUTINES 

The  original  version  of  LARM2  contained  heat  as  the  only 
constituent  transported  because  of  the  dependence  of  the  dynamics 
on  density.  Heat  was  transported  with  the  laterally  averaged 
advective-dispersive  relationship  in  which  attenuated  short-wave 
solar  radiation  and  surface  heat  exchange  were  included  as 
source-sink  terms.  The  transport  relationship  was  solved  im¬ 
plicitly  layer  by  layer  using  the  computationally  efficient 
Thomas  algorithm. 

Use  of  LARM2  for  the  study  of  the  transport  of  an  additional 
constituent,  residual  chlorine,  (Edinger  and  Buchak,  1978),  was 
achieved  by  using  the  advective-dispersive  transport  relationship 
with  residual  chlorine  source  and  sink  terms.  The  computations 
were  performed  by  replicating  the  advective  and  dispersive  terms 
of  the  heat  transport  relations  and  writing  the  appropriate 
sources  and  sinks.  The  LARM2  simulations  thus  solved  the  trans¬ 
port  equations  twice  on  each  time  step:  once  for  heat,  and  a 
second  time  for  residual  chlorine.  Adding  a  complete  set  of 
transport  computations  for  the  second  constituent  almost  doubled 
the  computational  time  of  a  LARM2  simulation. 

Studies  of  the  LARM2  code  showed  that  about  30%  of  the  sim¬ 
ulation  time  was  devoted  to  the  transport  relationships, and  most 
of  this  time  was  utilized  in  evaluating  the  tridiagonal  transport 
coefficients.  The  computations  could  be  made  more  efficient  for 


additional  water  quality  constituents  if  the  transport  coeffi¬ 
cients  could  be  evaluated  only  once  per  iteration  and  used  re¬ 
peatedly  for  additional  constituents. 

One  method  studied  to  generalize  the  transport  equations 
for  additional  constituents  was  to  write  them  in  matrix  form 
and  invert  the  matrix  at  each  time  step.  The  inverse  matrix  then 
multiplies  the  constituent  source-sink  terms  to  give  constituent 
concentrations.  The  matrix  inversion  process  took  almost  ten 
times  the  computational  time  for  tridiagonal  evaluation  of  the 
transport  equation  and  required  storage  that  increased  as  the 
square  of  the  number  of  active  cells.  Standard  matrix  inversion 
algorithms  also  did  not  retain  significant  accuracy  to  maintain 
proper  heat  and  mass  balances. 

The  second  and  adopted  method  was,  simply,  to  evaluate  the 
tridiagonal  coefficients  for  all  lines  of  the  reservoir  grid  and 
retain  them  for  the  successive  constituent  computations.  This 
method  of  computation  could  easily  be  generalized  to  a  water 
quality  transport  module  (WQTM)  that  carries  out  the  transport 
computations  for  each  water  quality  constituent  after  evaluation 
of  sources  and  sinks. 

Formalization  of  WQTM 

A  formal  statement  of  the  WQTM  can  be  derived  from  the 
laterally  averaged,  vertically  integrated  advective-dispersive 
transport  relationship  for  any  constituent,  C,  regardless  of  its 
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sources  or  sinks.  The  formalization  shows  that  the  transport 
coefficients  need  be  evaluated  only  once,  and  only  the  source- 
sink  terms  need  be  evaluated  for  each  constituent. 

The  laterally  averaged,  vertically  integrated  constituent 
transport  relationship  is 


a!  «hc>  +  sf  ‘"Bhc>  +  <»bbc)w%  '  <W 


.(,M,  +  (D  2S£, 

9x  '  3z  ;k+»s  **  3z  ;k-Ji 


Hn  Bh 
V 


(3.1) 


where 

b  laterally  averaged  lake  width  at  top  or  bottom  of  cell 
face  (m) 

B  laterally  averaged  lake  width  integrated  over  h  (m) 

C  laterally  averaged  constituent  concentration  integrated 
over  h  (mg  n-1 ) 

D  x-direction  temperature  and  constituent  dispersion 
coefficient  (m2/s) 

D  z-direction  temperature  and  constituent  dispersion 

z  coefficient  (m2/s)‘ 

h  horizontal  layer  thickness  (m) 

Hq  source  strength  (mg- i-1.  m3»  s-1  ) 

k  integer  layer  number,  positive  downward 

t  t ime  ( s ) 

U  x-direction,  laterally  averaged  velocity  integrated 
over  h  (m/s) 

V  cell  volume  (B*h-Ax)  (m3). 

wh  z-direction,  laterally  averaged  velocity  (m/s) 

D  (W  in  FORTRAN  code) 

x  and  z  Cartesian  coordinates:  x  is  along  the  lake  center- 
line  at  the  water  surface,  positive  to  the  right,  and  z 
is  positive  downward  from  the  x-axis  (m) 
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In  LARM2  the  velocities,  U  and  W,  and  dispersion  coeffi¬ 
cients,  Dv  and  D_,  are  evaluated  from  the  hydrodynamic  equations 
using  Richardson  number  relationships  and  are  available  for 
evaluating  the  constituent  balances.  Each  constituent  being 
evaluated  has  a  balance  given  by  Equation  3.1,  and  each  balance 
must  be  evaluated  separately.  Reactions  and  interactions  be¬ 
tween  constituents  are  included  in  the  source-sink  term,  Hq, 
which  is  evaluated  separately  from  the  transport  equations. 

For  laterally  averaged  vertically  integrated  transport 
where  the  horizontal  grid  length  is  much  greater  than  the  verti¬ 
cal  grid  length,  the  constituent  transport  relationships  can  be 
evaluated  layer  by  layer  along  the  x-axis  in  an  implicit  finite 
difference  formulation  on  an  i,k  grid  with  i  horizontal  and  k 
vertical.  The  upwind  differencing  scheme  used  in  the  heat  bal¬ 
ance  computation  to  maintain  exact  conservation  without  averaging 
concentrations  is  also  applied  to  the  constituent  balance. 

The  implicit  tridiagonal  form  of  the  transport  equation  is 

maintained  by  taking  the  longitudinal  transport  terms  forward  in 

time  and  lagging  the  vertical  transport  terms.  The  relationship 

could  be  expanded  to  an  alternating-direction- implicit  scheme 

(ADI)  by  evaluating  the  left-hand  side  for  the  x-direction  terms 

over  half  a  time  step  and  then  evaluating  the  vertical  terms  on 

the  left-hand  side  for  the  second  half  a  time  step.  The  ADI  is 

unnecessary  in  LARM2  because  of  the  much  greater  size  of  Ax 

compared  to  h,  and  it  is  incompatible  since  D  is  known  only  at  a 

z 


lagged  time  step  and  must  be  evaluated  with  the  proper  gradient, 
3C/3z,  at  that  time  step.  A  third  reason  for  lagging  the  ver¬ 


tical  transport  terms  is  that  Dz  may  be  different  from  one  con¬ 
stituent  to  the  next,  and  it  is  easily  evaluated  as  part  of  the 
right-hand  side  of  the  finite  difference  form  of  Equation  3.1. 

The  spatially  implicit  transport  relationship  can  be  ex¬ 
pressed  in  tridiagonal  form  as : 


a.C!  i  l.  +  v.c!  ,  +  c,  Cjj.1  .  =  d. 
i  i-l,k  i  i,k  i  i+l,k  i 


for  each  k  layer  where: 
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where 

Aik*  right  cell  face  area  =  [(Bh)i+1  k 
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US^  k  =  indicator  of  flow  direction  for  upstream  differencing 

*  1  if  Ui  k  >  0,  otherwise  US^  k  *  0 

USi-l,k  =  1  if  Ui-l,k  "  °’  otherwise  =  0 

Equation  3.2  shows  that  each  layer  on  the  LARM2  grid  has  a  set 
of  four  tridiagonal  coefficients,  a,  v,  c,  and  d,  on  each  time 
step.  Furthermore,  three  of  these,  a,  v,  and  c,  depend  only  on 
the  velocities,  geometry,  and  dispersion  coefficients  and  are 
invariant  in  a  time  step  from  constituent  to  constituent.  The 
transport  coefficients,  a,  v,  and  c,  need  be  evaluated  only  once 
per  time  step,  saved,  and  used  for  all  other  water  quality 
constituents.  In  LARM2,  the  transport  coefficient  arrays,  a,  v, 
and  c,  are  evaluated  where  they  are  used  first  in  the  heat  bal¬ 
ance  and  then  retained  for  use  with  the  remaining  water  quality 
constituents . 

The  tridiagonal  coefficient  d  is  dependent  on  the  partic¬ 
ular  constituent  being  evaluated.  It  includes  the  storage  term 
or  old  concentrations.  Secondly,  it  includes  the  two  vertical 
transport  terms,  advection  and  mixing,  which  are  lagged  in  time. 
Lastly,  it  includes  the  source-sink  term,  HQ,  which  must  be 
evaluated  for  each  constituent  reaction  and  interaction.  It  is 
convenient  to  design  WQTM  to  evaluate  the  d  coefficient  terms 
for  each  constituent  over  each  time  step  and  to  call  the  tri¬ 
diagonal  solver  (subroutine  TRIDAG)  for  that  constituent. 


The  FORTRAN  coding  of  the  WQTM  algorithm  is  shown  in  Table 
3.1.  Constituent  concentrations  at  the  new  time  step  for  con¬ 
stituent  JC  are  computed  on  each  pass  through  the  algorithm 
(DO  loop  1820).  The  algorithm  begins  with  the  evaluation  of  the 
source-sink  term,  Hq,  which  is  followed  by  the  retrieving  of  the 
a,  v,  and  c  vectors  from  their  storage  arrays  and  the  assembling 
of  the  d  vector  from  its  components,  including  Hn>  and  finally 
the  use  of  the  tridiagonal  solution  algorithm.  The  steps  follow¬ 
ing  the  evaluation  of  Hn  are  performed  for  each  layer  in  the  grid 
from  top  to  bottom.  Note  that  the  top  layer  d  computation  uses 
a  separate  set  of  vertical  velocities  that  are  computed  from  the 
change  in  mass  storage  in  the  top  layer,  rather  than  from  the  con¬ 
tinuity  expression  around  each  cell  in  the  next  lower  layer  which 
is  the  case  for  the  every  other  layer.  This  procedure  is  used  to 
maintain  perfect  constituent  balances  and  is  taken  from  the  heat 
balance  computation. 

The  evaluation  of  H  begins  with  the  initialization  of  the 

n 

HQ  array,  since  one  array  is  used  for  all  constituents  and  the 
heat  balance.  Secondly,  HQ  is  augmented  by  the  reaction- 
interaction  rates  for  the  current  constituent  in  terms  of  every 
other  constituent  and  all  other  internal  sources  and  sinks  (decay, 
settling,  etc.).  Finally,  external  sources  and  sinks  (inflows 
and  withdrawals)  are  considered. 

Constituent  Internal  Sources  and  Sinks 

The  internal  sources  and  sinks  evaluation  in  WQTM  includes 
source-sink  and  reaction  terms  for  each  biochemical  water  quality 
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parameter  being  evaluated.  Each  constituent  has  reaction  rates, 
settling  velocities,  etc. ,  and  interactions  with  other  constit¬ 
uents  that  are  evaluated  in  WQTM. 

The  evaluation  of  HQ  in  WQTM  is  based  on  the  fact  that  the 
transport  coefficients  are  independent  of  concentration  (after 
evaluation  for  density-dependent  terms)  and  that  almost  every 
constituent  reaction  can  be  written  in  the  form: 


6C1 

6t 


-  -K1C1+K2 


CIC2 

a+BCl 


+ 


(3.7) 


where  <5Cl/6t  represents  all  of  the  transport  and  storage  terms 
about  a  point  that  are  included  in  a,  v,  and  c.  The  Cl  and  C2 
are  concentrations  of  constituents  1  and  2.  The  remaining  terms 
are  rate  coefficients,  cycle  limits,  and  other  terms  representing 
reactions  and  interactions. 

After  the  constituent  reaction  relationships  are  developed 
by  the  user  as  above,  then  the  reaction  source-sink  term  becomes 
for  Cl: 


HN(I,K)  -  HN(1,K)  +  BH2(I,K) *DLX*RR(1 ,2)*C2(I,K,2)  +  ...  (3.8) 

where  RR(1,2)  is  the  rate  at  which  constituent  1  is  transferred 
to  constituent  2  and  C2(1,K,2)  is  the  concentration  of  constit¬ 
uent  2  at  the  old  time  level.  The  H  term  needs  to  be  evaluated 

n 

within  the  format  of  the  LARM2  geometry  computations,  and  this 
is  another  function  of  WQTM. 

The  user-specified  portion  of  WQTM  requires  knowing  the 


number  of  constituents  being  considered,  the  rate  expressions  as 
given  in  the  form  of  Equation  3.7,  and  how  the  rate  parameters 
are  evaluated.  It  is  the  last  that  is  written  into  the  user- 
defined  portion  of  WQTM  for  each  constituent.  The  procedure  for 
assembling  the  user  portions  of  the  WQTM  is  as  follows:  First 
write  the  rate  expressions  for  the  problem  being  considered  and 
expressions  for  the  controlling  parameters;  this  can  usually  be 
done  in  traditional  format  before  translating  to  WQTM  notation. 
Second,  complete  the  user-specified  statements.  The  procedure 
will  be  demonstrated  for  a  few  examples. 

Example  -  Sediment  Transport 

Consider  a  sediment  of  narrow  size  range  whose  settling 
rate  is  described  by  a  single  settling  velocity,  Vs.  The  sediment 
may  also  be  scoured  at  the  bottom  p.t  a  rate  proportional  to  the 
adjacent  horizontal  velocity. 

From  the  surface  to  the  bottom,  the  local  change  in  con¬ 
centration  is 

-  Vs  3C/3z  (3.9) 

The  effect  of  the  vertical  advective  velocity,  W(I,K),  is  al¬ 
ready  accounted  for  in  the  WQTM  advective  transport.  For  bottom 
scour  and  resuspension,  the  shear  function  SF  is  used,  and  for 
the  bottom  cells: 

3C 
at 


SF*U(I,KB)  -  Vs  3C/3z 


(3.10) 


The  reaction-interaction  sources  and  sinks  become  for  the  sur¬ 
face  layer: 

HN(I,KT)  =  HN(I,KT)  -  BH2(I,KT)*DLX*[Vs*C2(I,KT, 1) ]/H2(I,KT)  (3.11) 

for  the  bottom  layer: 

HN(I,KB)  =  HN(I,KB)  +  BH2(I,K)*DLX*[SF*U(I,KB) ] 

-  Vs*[C2(I,KB,l)  -  C2(I,KB-1,1)/H2(I,KB)]  (3.12) 

and  for  internal  layers: 

HN(I,K)  “  HN(I,K) 

-  BH2(I,K)*DLX*Vs*[C2(I,K, 1)  -  C2(I,K-1,1) ]/H2(I,K)  (3.13) 

These  are  inserted  in  the  WQTM  routine  for  this  constituent. 

WQTM  then  evaluates  the  transport  of  the  constituent  and  gives 
the  solution  vector.  The  sediment  inflow  concentrations  are 
specified  as  data  input  for  evaluation  of  the  transport  source 
and  sink  contributions  in  WQTM. 

The  sediment  transport  internal  source-sink  routine  is 
shown  for  a  single  "sediment  concentration"  as  an  example  only. 

With  the  efficient  transport  computations  as  provided  in  WQTM, 
it  is  more  realistic  to  compute  the  transport  of  a  number  of 
ranges  of  sediment  sizes,  each  range  having  its  specific  settling 
velocity  and  bottom  scour  functions,  and  possibly  interactions 
between  size  ranges.  LARM2  is  sufficiently  flexible  to  include 
the  effects  of  sediment  concentrations  on  density  and,  subsequently, 
on  the  velocity  field. 


The  constituent  internal  sources  and  sinks  can  be  used  to 
set  up  the  interacting  rate  expressions  for  any  number  of  con¬ 
stituents,  once  those  expressions  are  known.  One  of  the  more 
commonly  understood  multiconstituent  water  quality  processes  is 
the  nitrogen  cycle  in  nitrogen-limiting  systems.  Nitrogen 
limiting  means  that  there  is  sufficient  phosphorus  available  so 
as  not  to  control  plant  growth  and  so  as  to  allow  all  stages  of 
the  nitrogen  cycle  to  develop  over  a  season. 

One  description  of  a  seven-stage  nitrogen  cycle  has  been  de¬ 
veloped  by  Najarian  and  Harleman  (1975).  The  cycle  is  shown  in 
Figure  3.1  and  has  seven  compartments  of  nitrogen  (N),  including 
(1)  ammonia;  (2)  nitrite;  (3)  nitrate;  (4)  phytoplankton  nitro¬ 
gen;  (5)  zooplankton  nitrogen;  (6)  particulate  organic  nitrogen; 
and,  (7)  dissolved  organic  nitrogen.  The  seven  rate  expressions 
for  each  compartment  are  found  from  the  paths  in  Figure  3.1. 

They  are: 
ammonia: 


"  R41N4  +  R51N5  +  R71N7 


N1N4 

R12N1  “  R14  K1+N1 
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nitrite : 
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nitrate : 
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zooplankton  nitrogen: 
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particulate  organic  nitrogen: 
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6T'R46N4  +R56N5-R67N6 


dissolved  organic  nitrogen: 


SN? 

Tt~  =  R47N4  +  R67N6  *  R71N7 


(3.16) 
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The  transfer  rates  R^j  are  complex  functions  of  other  variables, 
including  sunlight,  biomass,  temperature,  and  even  the  concen¬ 
tration  of  to  Ng.  The  transfer  rates  should  be  evaluated  in 
WQTM  ahead  of  the  source-sink  evaluation. 


The  growth-limiting  uptake  and  production  rates,  such  as 
the  reduction  of  ammonia  by  phytoplankton  (r14  N1N4/(K1  + 

in  Equation  3.14  and  its  transfer  to  phytoplankton  nitrogen  in 
Equation  3.17,  make  the  source-sink  terms  relatively  complex 
statements.  For  short  computational  time  steps,  these  processes 
can  be  linearized  into  the  rate  coefficients,  so  that  the  evalua¬ 
tion  of  HN ( I , K )  for  each  constituent  reduces  to: 

DO  JC  =  1,  NC 

HN(I,K)  =  HN(I,K)  +  RR( JC ,M) *C2 ( I , K , M) *BH2 ( I , K) *DLX 

CONTINUE 

The  overall  rate  coefficient  RR(JC,M),  multiplying  constituent 
M  to  get  constituent  JC,  is  summarized  in  Table  3.2  for  each 
relationship,  and  they  are  evaluated  prior  to  HN(I,K).  The  seven 
rate  expressions.  Equations  3.14  to  3.20,  are  reduced  to  an  easily 
evaluated  form  in  WQTM  in  terms  of  twenty  rate  expressions. 

An  alternative  scheme  is  to  write  out  each  relationship  of 
Equations  3.14  to  3.20  in  WQTM  for  evaluation  of  the  source-sink 
terms,  which  leads  to  a  more  lengthy  expression.  Since  the 
source-sink  terms  are  evaluated  from  concentrations  at  the  pre¬ 
vious  time  step,  the  linearized  summation  is  equivalent. 

The  nitrogen  cycle  can  be  coupled  to  a  dissolved  oxygen 
balance  that  has  a  sink  utilization  by  oxidation  of  ammonia  to 
nitrate,  as  well  as  production  and  respiration  by  plankton. 

Writing  the  internal  source-sink  routine  is  the  same  for  the 
D.O.  balance  as  for  the  nitrogen  balance:  (1)  begin  with  the 


basic  rate  expressions;  (2)  translate  to  a  source-sink  expression 
for  each  constituent.  Inflows  and  outflows,  as  well  as  transport 
through  the  waterbody,  are  evaluated  in  the  routine  WQTM  following 
the  source-sink  specifications. 
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Nitrogen-Cycle  Structure  in  Aerobic  Aquatic 
Ecosystems  (From  Najarian  and  Harleman  (1975)) 
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Summary  of  RR(JC,U)  Coefficients  for 
Linear  Evaluation  of  Nitrogen  Cycle 
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4.  SPECIAL  TOPICS  IN  LARM2  DEVELOPMENT 


The  examination  and  testing  of  numerous  longitudinal- 
vertical  hydrodynamics  codes  by  Johnson  (1981),  including 
earlier  versions  of  LARM2,  have  indicated  a  number  of  topics 
that  require  further  consideration.  These  include  other  pos¬ 
sible  geometric  configurations,  refinements  of  the  hydraulic 
computations,  and  other  methods  of  computing  turbulent  mixing 
parameters . 

The  topics  examined  in  the  context  of  the  present  LARM2 
development  are:  (1)  coordinate  transforms  and  bottom  slopes; 
(2)  variable  longitudinal  grids;  (3)  steady-state  solutions; 

(4)  channel  conveyance;  (5)  reservoir  branching;  and  (6)  tur¬ 
bulence  and  mixing  processes.  The  structure  of  the  LARM2 
theory,  computational  algorithms,  coding,  and  development  were 
examined  to  determine  how  the  topics  applied  to  LARM2  and  how 
they  might  be  accommodated. 

Coordinate  Transforms  and  Bottom  Slopes 

The  WES  GRH  flume  experiments  with  LARM2  at  a  very  small 
grid  size  (Johnson,  1981)  showed  that  a  cold  water  density  under 
flbw  in  the  sloping  flume  moved  faster  than  computed  by  LARM2. 
Sensitivity  analyses  with  LARM2  showed  the  computed  density 
underflow  speed  to  be  relatively  insensitive  to  bottom  friction 
but  highly  dependent  on  initial  inlet  conditions,  as  are  most 
density  flow  problems. 


Another  possible  effect  on  the  difference  between  the  ob¬ 
served  and  computed  density  front  speed  was  thought  to  be  rep¬ 
resentation  of  the  flume  bottom  slope  by  the  step-wise  vertical 
grid  of  LARM.  Two  methods  suggested  for  a  more  explicit  rep¬ 
resentation  of  bottom  slope  in  longitudinal-vertical  dynamics 
are  (1)  using  a  coordinate  system  transformed  from  the  LARM2  rect¬ 
angular  grid,  or  (2)  including  bottom  slope  explicitly  in  the 
LARM2  surface  wave  and  momentum  relations  rather  than  implicitly 
computing  it  from  the  grid  configuration. 

A  transformed  grid  system  that  is  approximately  parallel 
to  the  bottom  slope  is  shown  in  Figure  4.1.  The  transformed 
z  coordinate  increases  with  distance  down  the  reservoir.  The 
transformed  coordinate  is  mapped  onto  the  rectangular  compu¬ 
tational  grid  by  a  series  of  gradient  relationships  that  would 
require  rewriting  most  of  the  basic  computations  in  LARM2. 

Mapping  of  waterbodies  by  irregular  grids  with  transformation 
to  a  rectangular  grid  has  been  very  successful  in  two-dimensional, 
vertically-mixed  hydrodynamic  problems.  In  longitudinal-vertical 
dynamics,  it  is  necessary  to  consider  the  vertical  variation  of 
pressure  and  horizontal  density  gradient  in  the  transformation. 
Fixed-coordinate  transformations  for  longitudinal-vertical  dy¬ 
namics  have  not  yet  been  tested, even  for  the  most  elementary 
cases. 

The  type  of  transformation  shown  in  Figure  4.1  that  would 
apply  to  longitudinal-vertical  dynamics  has  a  number  of 


limitations.  The  computations  are  providing  extensive  vertical 
detail  at  the  upper  end  of  the  reservoir  where  it  is  not  usually 
needed,  and  very  limited  vertical  detail  at  the  dam  where  it  is 
usually  needed.  The  transformed  grid  also  provides  complica¬ 
tions  in  adding  or  subtracting  grid  detail  as  the  free  water 
surface  rises  and  falls.  In  order  to  handle  this  problem,  it 
may  be  necessary  to  use  a  time-varying  Lagrangian  grid  trans¬ 
formation  which  is  significantly  more  complex  than  the  dynamic 
computations  on  a  rectangular  grid. 

Bottom  slope  is  determined  in  the  LARM  computations  directly 
from  grid  geometry.  Often,  in  fitting  a  given  Ax-by-h  grid  to  a 
waterbody,  the  actual  bottom  slope  is  under- or  overstated  at  the 
point  of  the  horizontal  velocity  computation.  The  manner  in 
which  bottom  slope  is  implicitly  included  and  how  it  can  be 
explicitly  stated  can  be  shown  using  the  primitive  LARM2  rela¬ 
tionships  of  horizontal  momentum: 


du  1  3PX„ 

9t  p  9x 


(4.1) 


where  F  is  the  sum  of  all  other  terms; 
vertical  momentum  (hydrostatic  approximation) 


3P  „ 
31“  p* 


(4.2) 


and  vertical  integrated  continuity: 


where  n  is  the  water  surface  elevation  and  G  is  the  reservoir 


bottom  elevation.  Utilizing  the  vertical  pressure  distribution, 
the  horizontal  pressure  gradient  becomes: 

fz 


1_  3p  _  3n  £  3p_  dz 

p  5x  -  ®9x  P  3x 


(4.4) 


'n 


and  is  described  in  terms  of  the  surface  elevation  as  presently 
used  in  LARM2. 

Direct  inclusion  of  bottom  slope  in  the  horizontal  momentum 
is  achieved  simply  by  writing: 


n  -  H  +  c 


(4.5) 


where  H  is  the  total  water  depth.  The  horizontal  pressure 
gradient  becomes: 


1  3P  3H  3g  £  3p_  dz 

"p3ji:"85x  85x“P  n3x 


(4.6) 


where  3g/3x  is  the  bottom  slope.  So  long  as  3g/3x  is  determined 
from  the  grid  geometry,  then  the  horizontal  pressure  gradient 
by  Equation  4.6  is  identical  to  the  LARM2  form  in  Equation  4.4. 
However,  Equation  4.6  states  that  the  bottom  slope  can  be  eval¬ 
uated  independently  of  the  grid.  It  is  then  only  necessary  to 
rewrite  the  horizontal  pressure  gradient  evaluation  in  terms  of 


tbe  gradient  of  water  column  depth  rather  than  water  surface 
slope.  The  surface  equation  readily  translates  to  a  depth  com¬ 
putation  using  3n/3t  =  3H/ 3t,  and  bottom  slope  is  included  as 
one  of  the  forces. 

Variable  Longitudinal  Grid 

The  finite  difference  formulations  in  LARM2  are  presently 
developed  for  a  uniform  longitudinal  spacing  (constant  Ax). 
Usually,  in  reservoir  and  estuary  problems,  the  grid  spacing  con¬ 
veniently  is  of  the  order  of  1  km  to  2  km.  In  reservoir  prob¬ 
lems  such  as  withdrawal  and  pumpbacks  or  locations  of  complex 
tributary  and  geometry  configurations,  it  is  sometimes  useful 
to  have  more  spatial  detail  of  the  reservoir  flow  field  and 
transport . 

Any  finite  difference  scheme  of  the  momentum  and  transport 
equations  can  be  examined  for  consistency  by  identifying  the 
"cell"  around  the  primary  variable  being  computed.  The  momentum 
balance  can  be  examined  for  a  cell  around  the  location  of  the 
horizontal  velocity  component,  and  the  transport  balance  can  be 
examined  for  a  cell  around  the  location  of  the  constituent 
concentration.  Each  cell  has  gradient  at  both  faces,  which  are 
the  surface  slope  and  horizontal  density  gradients  for  momentum 
and  the  advection  and  dispersion  gradients  for  constituent 
transport . 

The  constituent  transport  relation  enters  the  momentum  bal¬ 


ance  through  the  horizontal  density  gradient.  For  a  space- 


staggered  grid  where  constituent  concentrations  are  computed  at 
points  between  the  velocities,  the  gradients  of  velocity  in 
transport  are  compatible  to  the  gradients  of  density  in  momen¬ 
tum  when  the  grid  spacing  is  uniform.  The  same  is  true  for  the 
water  surface  profile  equation. 

Unequal  grid  spacing  requires  extensive  averaging  of  pri¬ 
mary  variables  by  weighting  between  grid  points.  The  conven¬ 
ient  upwind  differencing  on  momentum  and  constituent  advection 
inherently  assumes  a  uniform  grid  spacing  and  would  require  a 
weighted  form  for  a  nonuniform  grid  to  be  compatible  with 
continuity.  It  is  possible  to  develop  the  finite  difference 
forms  of  the  relationships  for  an  unequal  grid  spacing,  but 
they  must  obey  the  same  conditions  of  consistency,  compatibility, 
and  continuity  as  found  for  the  uniform  spacing. 

The  perceived  limitations  of  the  uniform  Ax  used  in  LARM2 
can  be  overcome  by  two  methods.  First  is  to  perform  computa¬ 
tions  at  a  Ax  smaller  than  normally  used.  The  LARU2  computa¬ 
tions  are  quite  efficient  in  terms  of  computer  time,  yet  most 
problems  are  run  with  15-20  longitudinal  segments  and  up  to 
25  vertical  layers.  The  number  of  longitudinal  segments  could 
easily  be  doubled  without  encountering  excessive  computer  costs. 

Another  method  is  to  perform  LARM2  simulations  for  the 
usual  large  Ax  and  then  set  up  a  second  LARM2  simulation  for  a 
portion  of  the  waterbody  at  a  smaller  Ax.  The  computed  veloc¬ 
ities,  concentrations,  etc.,  of  the  large-scale  case  become  the 


boundary  conditions  for  the  smaller  scale  case.  This  approach 
presumes  there  is  similarity  over  the  scales  of  different  Ax, 
which  there  is  for  a  uniform  grid  spacing. 

Steady-State  Solution 

There  are  a  few  problems  in  longitudinal-vertical  water- 
body  dynamics  that  can  be  examined  using  steady-state  solutions. 
Steady- state  cases  have  steady  boundary  flows,  steady  boundary 
consitutent  concentrations,  and  steady  boundary  exchange 
processes.  Steady-state  conditions  might  be  specified  for  pre¬ 
liminary  examination  of  flow  fields  before  a  complete  set  of 
time- varying  boundary  data  is  assembled.  The  real  reservoir 
problem  is  one  of  time-varying  unsteady  inflows,  outflows,  and 
meteorological  conditions. 

The  LARM2  computations  are  designed  to  iterate  over  time 
the  sequence  of  the  surface-wave  equation,  the  pressure  distri¬ 
bution,  the  horizontal  momentum  balance,  internal  continuity, 
tations  can  be  performed  for  specified  steady  boundary  condi¬ 
tions  which  are  really  a  special  case  of  the  more  general 
time-varying  boundary  conditions.  For  iterative  solutions 
of  the  time-varying  equations  to  steady  state,  the  flow  field 
establishes  rather  quickly,  while  the  constituent  transport  takes 
considerably  longer  to  establish  the  constituent  distribution. 

In  general,  for  an  implicit  solution,  the  flow  field  is  estab¬ 
lished  within  two  surface-wave  travel  times  over  the  length  of 


the  reservoir  or  estuary,  while  the  constituent  balance  gets 
established  over  the  residence  time  of  the  waterbody.  The  latter 
is  due  to  the  time  required  to  build  up  storage  of  the  constit¬ 
uent  within  the  waterbody. 

Steady- state  solutions  appear  attractive  because  of  the  re¬ 
duction  in  computer  time.  However,  even  if  the  local  time  change 
components  are  eliminated  from  the  formulation  of  the  basic 
equations,  the  solution  technique  must  be  iterative.  For  steady- 
state  conditions,  the  LARM2  computations  can  be  made  quite  effi¬ 
cient  by  eliminating  the  local  storage  term  from  only  the  con¬ 
stituent  transport  relation.  It  is  done  quite  simply  by  elimin¬ 
ating  the  BH1( I , K)/DLT  from  the  V(I)  of  the  tridiagonal 
coefficient  and  T2(I ,K)*BH2(I ,K)/DLT  term  from  the  D(I)  tri¬ 
diagonal  coefficient  in  the  constituent  transport  equations. 

These  changes  produce  the  steady- state  constituent  balance;  the 
computational  procedures  remain  unchanged. 

A  LARM2  steady-state  solution  is  achieved  basically  by 
computing  a  steady-state  constituent  distribution  for  each 
iteration  of  the  flow-field  dynamics.  It  is  based  on  the  fact 
that  the  implicitly  computed  flow  field  becomes  established 
quite  rapidly  and  without  initialization  oscillations.  The 
computational  procedure  is  similar  to  iterative  solution  tech¬ 
niques  that  would  be  required  by  reformulation  from  the  steady- 
state  form  of  the  basic  equations. 


Channel  Conveyance 


Channel  conveyance  refers  to  the  fact  that  a  larger  fraction 
of  flow  per  cross-sectional  area  (or  velocity)  occurs  in  main 
channel  areas  than  occurs  in  overbank  areas.  For  LARM2  this 
means  choosing  the  lateral  widths  to  apply  to  the  main  channel 
area  and  treating  the  overbank  areas  as  regions  of  lateral  in¬ 
flow  and  outflow.  The  result  is  generally  higher  velocities  in 
the  main  channel  area  than  would  be  obtained  using  widths  for 
the  whole  cross  section.  Methods  of  computing  the  channel  con¬ 
veyance  sections  and  widths  are  presently  available  in  the 
Hydrologic  Engineering  Center's  GEDA  program,  once  the  user  has 
specified  which  channel  regions  are  to  be  included  in  the  con¬ 
veyance  area. 

Revising  LARM2  to  incorporate  a  conveyance  width  and  an 
overbank  area  can  be  performed  within  the  existing  computational 
structure.  At  each  longitudinal  location,  an  overbank  planar 
area,  AB(I,K),  can  be  introduced  which  is  a  function  of  depth. 
Presently,  a  tributary  inflow,  QTRIB(J),  is  specified  for  each 
J  tributary.  Each  J  tributary  has  longitudinal  position,  ITRIB, 
and  a  vertical  position,  KTRIB,  computed  on  the  basis  of  density 
inflow  depth.  The  QTRIB(J)  is  presently  the  lateral  inflow  to 
the  mainstem  LARM2  computations.  It  can  be  used  to  account  for 
lateral  inflows  and  outflows  to  and  from  the  main  flow  as: 


l 


QTRIB(J)  -  QTRIB(J)  -  AB(I,K)  |^ 


(4.7) 


where  3n/3t  is  the  change  in  water  surface  elevation  from  one 
time  step  to  the  next.  The  overbank  flow  computation  indicated 
by  Equation  4.7  assumes  that  the  vertical  distribution  of  tem¬ 
perature  and  other  water  quality  parameters  is  the  same  in  the 
tributary  overbank  area  as  in  the  mainstem  area.  It  also  assumes 
that  the  net  flow  into  and  out  of  the  overbank  area  when  there 
is  no  tributary  inflow  is  due  only  to  change  in  water  level. 

The  above  correction  allows  using  a  more  realistic  conveyance 
width  in  the  main  LARH2  computations  past  mouths  of  tributary 
embayments.  Larger  tributary  segments  can  be  handled  by 
branching. 

Reservoir  Branching 

The  LARM  code  is  presently  structured  to  compute  vertical 
velocity  and  constituent  profiles  along  a  single  center  line, 
with  tributary  inflows  and  withdrawals.  The  case  of  a  reservoir 
formed  by  a  dam  near  the  junction  of  two  major  tributaries  has 
been  handled  successfully  in  LARM2  by  running  the  model  center 
line  continuously  down  one  arm  and  up  the  other.  More  dendritic 
reservoir  geometries  have  a  number  of  major  branches  for  which 
it  is  necessary  to  have  longitudinal  and  vertical  resolution  of 
the  velocity  and  constituent  fields. 

The  LARM2  code  is  presently  formulated  for  "flow-flow" 
boundaries  in  which  inflows  and  outflows  at  either  end  of  the 
model  are  specified.  The  branching  problem  with  a  major  branch 


4.11 


intersecting  the  mainstem  reservoir  is  an  "inflow-head"  boundary 
problem, with  inflows  specified  in  the  upper  end  of  the  branch 
and  the  water  surface  matched  at  the  junction.  LARM2  has  been 
run  successfully  for  the  "inflow-head"  boundary  case  in  the  es¬ 
tuary  version,  LAEM,  (Edinger  and  Buchak,  1980).  The  proper 
form  of  the  dynamic  boundary  conditions  and  their  location  in 
the  computation  are  known  and  have  been  demonstrated. 

Extension  of  LARM2  to  computing  dynamics  in  branches  can 
be  accomplished  by  computing  on  each  time  step  the  mainstem 
dynamics  with  flow-flow  boundaries,  then  computing  the  flow- 
head  dynamics  for  each  branch  using  the  mainstem  elevation.  The 
inflows  and  outflows  between  the  tributary  branch  and  mainstem 
are  computed  for  the  branch  dynamic  computations.  These  flows 
become  tributary  flows  to  the  main  branch. 

The  LARM  computational  algorithms  need  few  changes  to  in¬ 
corporate  branching  cases.  The  longitudinal  computational  limits 
have  been  generalized  to  a  variable  beginning  and  ending 
coordinate.  The  end  coordinates  can  be  made  a  function  of  the 
branch  numbei .  For  a  typical  case  the  mainstem  may  run  from 
1=1  to  1=16,  the  first  branch  1=17  to  1=28,  and  the  second  branch 
1=29  to  1=34.  This  procedure  allows  using  the  present  arrays  of 
variables  and  computational  loops.  An  additional  coordinate  is 
required  to  specify  where  the  branches  intersect  the  mainstem 
and  where  to  apply  the  branch  boundary  head  condition.  The  pre¬ 
sent  computational  structure  of  LARM2  is  such  that  extension  to 
branching  problems  is  quite  feasible. 


Turbulence  and  Mixing  Processes 

The  computational  algorithms  in  LARM2  are  developed  to  in¬ 
clude  turbulent  transport  dispersion  coefficients  and  eddy  vis¬ 
cosities  as  functions  of  time  and  space.  The  four  parameters 
that  are  utilized  for  turbulent  transport  of  constituent  and 

momentum  are  the  vertical  dispersion  coefficient,  D  ,  the  ver- 

z 

tical  eddy  viscosity,  A,,  tne  longitudinal  dispersion  coeffi- 

z 

cient,  D  ,  and  the  longitudinal  eddy  viscosity,  A  . 

LARM2  utilizes  the  Richardson  number  concept  to  evaluate 
the  vertical  dispersion  coefficient  and  vertical  eddy  viscosity 
as  functions  of  buoyancy  and  velocity  shear.  It  is  a  classical 
formulation  disscussed  in  Edinger  and  Buchak  (1980)  and  has  the 
form  of : 

Az  -  Azo  (1  +  10  Ri)”  1/2  (4.8) 


for  momentum  and: 


D-D  (1+10  Ri)"  3/2 
2  zo  — 


(4.9) 


for  constituent  transport.  The  Richardson  number,  Ri,  is 
defined  as: 


Ri 
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and  is  the  ratio  of  potential  energy  due  to  buoyancy  to  kinetic 


energy  being  dissipated.  In  the  classical  case,  AZo  and  Dzo 
are  taken  as  constants  for  the  waterbody.  The  Az  and  Dz  vary 
spatially  and  temporally  throughout  the  waterbody  and  are  con¬ 
strained  between  their  molecular  values  (as  Ri  gets  large)  and  a 
maximum  of  h2/At  for  the  grid  size  scale  effects.  The  longi¬ 
tudinal  coefficients  D  and  A  are  presently  taken  as  constants 

X  X 

since  the  solutions  are  insensitive  to  them  at  large  scales. 

The  Richardson  number  formulation  accounts  for  changes  in 
turbulent  dispersion  and  eddy  viscosity  under  stratified  con¬ 
ditions  and  is  quite  simple  to  apply.  It  does  not,  however, 

allow  for  varying  A  and  D  for  unstratified  conditions  or  for 

z  z 

any  transport  of  turbulent  kinetic  energy  from  one  portion  of 

the  waterbody  to  another.  Nor  does  it  allow  for  similarity 

relations  between  D  and  D_  or  A  and  A  in  terms  of  modelling 

X  Z  X  z 

scales  Ax  and  Az,  except  empirically. 

Another  method  for  relating  the  dispersion  and  viscosity 
parameters  to  the  mean  flow  field  is  through  the  evaluation 
and  transport  of  turbulent  kinetic  energy,  K,  as  generated  by 
shear  and  buoyancy.  The  turbulent  kinetic  energy  is  a  scalar 
quantity  which  is  transported  as  a  constituent.  For  laterally 
averaged  dynamics,  the  turbulent  kinetic  energy  transport 
relationship  is  (Rodi,  1980): 


3(BAz3K/3z 


BAz(3U/3z)2  +  gB  D2^|^  -  E 


(4.11) 


where  B  is  the  lateral  width.  The  left-hand  terms  are:  the 
local  change  in  turbulent  kinetic  energy,  advection  with  the 
mean  flow  field,  and  turbulent  transport.  The  right-hand  side 
is  the  production  of  turbulent  kinetic  energy  by  mean  velocity 
shear  and  by  buoyancy  and  its  dissipation,  E,  by  viscosity. 

The  turbulent  dispersion  and  eddy  coefficients  are  related 
to  the  turbulent  kinetic  energy  as  (Rodi,  1980): 


A  -  C  L  K 
x  xx 


A  -  CL  K 
z  z  z 


(4.12) 


(4.13) 


where  C  and  C  are  constants  and  and  L  are  scale  lengths 
x  z  x  z 

related  to  the  size  of  the  waterbody.  It  is  thought,  (Edinger 

and  Buchak,  1980),  that  L  and  L_  are  related  to  the  size  of 

X  z 

the  computational  grid  in  numerical  models.  A  similar  set  of 

relationships  holds  for  and  D  . 

x  z 

Turbulent  kinetic  energy  dissipation,  E,  is  also  related 


to  K  as: 


E  -  C£  K3/2/L 


(4.14) 


where  Cj,  is  a  constant  and  L  is  a  scale  length.  It  could  also 
be  transported  similarly  to  K,  resulting  in  smaller  scales  of 
dissipation. 


If  turbulent  kinetic  energy  is  not  transported  or  is  taken 
as  zero,  the  right  hand  side  of  Equation  4.11  reduces  to  a 
Richardson  number  description  of  dissipation  similar  to 
Equation  4.8.  The  Richardson  number  formulations  thus  apply 
to  steady  flows  with  no  longitudinal  or  vertical  velocity 
variations.  This  is  seldom  the  case  in  reservoir  problems  and 
never  the  case  in  tidal  estuaries.  The  dispersion  coefficient 
and  eddy  viscosity  respond  immediately  to  buoyancy  and  shear  in 
the  Richardson  number  formulation,  while  evaluation  from  turbulent 
kinetic  energy  results  in  time  delays  between  the  occurrence  of 
shear  and  buoyancy  and  dispersion. 

Introducing  the  computation  of  turbu'ent  kinetic  energy  as 
generated  by  shear  and  buoyancy  and  dissipated  by  viscosity  into 
the  LARM2  code  is  a  relatively  easy  task  since  the  transport  com¬ 
putations  have  been  generalized  in  WQTM.  Its  use  allows  inves¬ 
tigating  higher-order  turbulence  transport  relationships  and 
examining  the  relationship  between  numerical  computational 
grid  scales  and  turbulent  length  scales. 

Another  value  of  the  turbulent  kinetic  energy  formulation  is 
in  evaluation  of  turbulence  in  surface  layers  due  to  wind  shear. 
Use  of  turbulent  kinetic  energy  transport  allows  wind-generated 
shear  to  be  transported  through  the  water  column  in  a  less  em¬ 
pirical  manner. 

The  turbulent  kinetic  energy  transport,  K,  is  computed  on 
the  scale  of  Ax,  Az,  and  B  for  a  given  location.  Reservoir 


withdrawals  and  pumpbacks  may  result  in  withdrawal  zones  or 
jets  which  are  at  a  smaller  scale  and  are  producing  or  dissi¬ 
pating  turbulent  kinetic  energy.  These  become  sources  or  sinks 
of  K  which  are  incorporated  in  its  transport  computation.  The 
turbulent  kinetic  energy  transport  formulation  allows  a  direct 
method  for  incorporating  turbulence  caused  by  withdrawal  and 
pumpback  directly  into  reservoir  analyses. 


Possibl 
of  Long! 


5 .  RECOMMEND AT IONS 


Recommendations  for  extension  of  LARM2  to  additional  areas 
of  study  are:  (1)  inclusion  of  channel  conveyance;  (2)  devel¬ 
opment  of  algorithms  for  reservoir  branching;  (3)  incorporation 
of  turbulent  kinetic  energy  transport  formulations  for  turbulent 
dispersion  and  eddy  viscosity  parameters;  and  (4)  summarization 
of  water  quality  constituent  internal  rate  expressions  into 
LARM2-WQTM  format. 

\ 

Channel  Conveyance 

Channel  conveyance  can  be  an  important  factor  in  reservoirs 
and  estuaries  that  have  extensive  overbank  areas  and  tributary 
embayments.  Incorporation  in  LARM2  requires  evaluation  of 
channel  conveyance  widths  and  modification  of  the  tributary 
flow  routines. 

Evaluation  of  the  channel  conveyance  widths  requires 
examination  and  possible  revision  of  the  EEC  GEDA-  program  to 
provide  conveyance  geometry  and  possibly  to  compute  overbank 
areas  as  a  function  of  elevation. 

Extension  of  the  LARM2  program  requires  (1)  modification  of 
the  tributary  inflow  routines  to  account  for  changes  in  overbank 
storage  and  (2)  modification  of  volume-area-elevation  computa¬ 
tions  to  account  for  overbank  volumes  and  areas.  All  of  these  re 
visions  can  be  made  with  the  present  structure  of  the  LARM2  code . 


Reservoir  Branching 

Reservoir  Branching  can  be  incorporated  into  the  present 
LARM2  code  by  the  addition  of  some  computational  algorithms.  The 
head-flow  type  of  boundary  condition  required  for  branching  has 
been  tested  in  the  estuarine  version  of  LARM2. 

Reservoir  branching  will  require  some  complex  coding  that 
will  require  careful  testing  as  it  is  developed  to  assure  con¬ 
tinuity  and  compatibility.  Also  required  is  revision  of  the 
LARM2  print  format  to  show  longitudinal-vertical  profiles 
in  each  arm. 

Branching  will  extend  the  capability  to  map  the  geometry 
of  reservoirs  onto  the  LARM2  computations  as  well  as  give  de¬ 
tailed  velocity  and  constituent  profiles  in  long  tributaries. 

Turbulent  Kinetic  Energy  Transport  and  Mixing 


The  vertical  and  longitudinal  dispersion  and  eddy  viscocity 
coefficients  can  be  evaluated  in  any  number  of  ways.  The 
present  formulations  in  LARM2  are  basically  steady-state 
Richardson  number  formulations.  The  generalized  transport  code 
in  LARM2  allows  using  the  transport  of  turbulent  kinetic  energy 
from  velocity  shear  and  buoyancy  as  a  basis  for  evaluating  the 
transport  processes. 

Turbulent  kinetic  energy  transport  is  being  used  for  eval¬ 
uation  of  dispersion  and  eddy  viscosity  coefficients  in  lon¬ 
gitudinal-vertical  estuarine  dynamics.  It  plays  the  important 
role  of  transporting  shear-generated  turbulence  away  from 
points  of  generation  where  artifically  high  eddy  viscocities 


might  be  computed  and  inducing  time  lags  in  the  velocity  field. 
Although  velocities  are  much  lower  in  reservoirs  than  in  es¬ 
tuaries,  the  vertical  velocity  shear  can  be  large,  partic¬ 
ularly  below  the  wind- driven  layers. 

The  LARM2  code  is  at  a  level  of  development  at  which  turbu¬ 
lent  kinetic  energy  transport  evaluation  of  dispersion  and  eddy 
viscosity  can  be  introduced  and  tested  without  uncertainity  of 
the  capacility  of  the  code  to  handle  it.  Its  further  development 
will  require  the  design  of  test  simulations  to  determine  the  ef¬ 
fects  of  including  turbulent  kinetic  energy  transport. 

Water  Quality  Consituent  Formulations 


The  WQTM  in  LARM2  has  been  structured  to  receive  the  inter¬ 
nal  source/sink  rate  expressions  for  any  number  of  interacting 
water  quailty  constituents.  The  latter  have  been  illustrated 
for  sediment  transport  with  settling  and  bottom  scour  and  for 
a  seven-constituent  nitrogen  cycle.  The  internal  source/sink 
computations  of  the  WQTM  are  general  and  can  be  written  for  any 
set  of  constituent  cycles  and  interactions  for  which  the  rate 
expressions  are  known. 

The  LARM2  notation  used  in  the  internal  source/sink  routine 
is  quite  simple;  and  the  structure  of  the  water  quality  expres¬ 
sions  in  the  routine  may  be  slightly  different  than  in  other 
computational  schemes,  such  as  the  one-dimensional  reservoir  of 


river  models.  It  is,  therefore,  recommended  that  the  water 


quality  and  quantitative  biology  routines  being  used  in  other 
Corps  programs  be  abstracted  and  summarized  in  LARM2/WQTM  format 
for  use  in  two-dimensional  reservoir  problems. 

There  are  numerous  water  quality  cycles  used  in  estuarine 
analyses  and  in  determining  the  fate  of  pollutants  that  are 
presently  not  among  the  Corps  water  quality  analysis  procedures. 
These  should  be  abstracted  from  the  literature  and  made  avail¬ 
able  in  LARM2/WQTM  format . 
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APPENDIX  A:  NOTATION 


Tridiagonal  matrix  coefficients 
Right  or  left  cell  face  area 

x-direction  momentum  dispersion  (i.e.  eddy  vis¬ 
cosity)  coefficient  (m2/s) 

z-direction  momentum  dispersion  (i.e.  eddy  vis¬ 
cosity)  coefficient  (m2/s) 

Neutral  stability  z-direction  momentum  disper¬ 
sion  coefficient  (m2/s) 

Laterally  averaged  lake  width  at  top  or  bottom 
of  cell  face  (m) 

Laterally  averaged  lake  width  integrated  over 
h  (m) 

B*h  (m2);  (BH1,  new  time  level,  BH2,  old  time 
level  in  FORTRAN  code) 

Chezy  resistance  coefficient,  m^/s 

Laterally  averaged  constituent  concentration 
integrated  over  h  (mg*fc-1)(Cl,  new  time 
level,  C2,  old  time  level  in  FORTRAN  code) 

Same  as  C  but  taken  at  the  new  time  level 

Resistance  coefficient 

Constant  in  turbulent  kinetic  energy  dissipa¬ 
tion  relation  (Equation  4.14) 

x-direction  temperature  and  constituent  disper¬ 
sion  coefficient  (m2/s) 

z-direction  temperature  and  constituent  disper¬ 
sion  coefficient  (m2/s) 

Neutral  stability  z-direction  temperature  and 
constituent  dispersion  coefficient 

Turbulent  kinetic  energy  dissipation 


f  * ' 


F 


Sum  of  other  terms  in  horizontal  momentum  equation 
(Equation  4.1) 

Acceleration  due  to  gravity  (m/s2) 


g 

h  Horizontal  layer  thickness  (m)(Hl,  new  time 

level,  H2,  old  time  level  in  FORTRAN  code) 

H  Total  depth  (m) 

Hq  Source  strength  for  heat  balance  (C*m3*s_1)  or 

constituent  balance  (mg • l~ 1 •m3s~ 1 ) (HN  in  FORTRAN 
code) 

i  Integer  segment  number,  positive  to  the  right 
(I  in  FORTRAN  code) 

j,JC,M  Index  to  denote  particular  water  quality  constituent 

J  Index  to  denote  particular  tributary 

k  Integer  layer  number,  positive  downward  (K  in 

FORTRAN  code) 

K  Turbulent  kinetic  energy 

Kl,K2,etc.  Water  quality  constituent  reaction  rates 

L  Scale  lengths  (related  to  waterbody  size) 

N  Nitrogen  stages 

P  Pressure  (Pa  =  N/m2) 

R  Rate  of  constituent  transfer,  s_l  (RR  in  FORTRAN 
code) 

Ri  Richardson  number 

SF  Shear  function  (mg* l~ 1 *m~ 1 ) 

t  Time  (s) 

T  Laterally  averaged  temperature  integrated  over 

h  (°C) 

u 


U 


x-direction  velocity  (m/s) 

x-direction,  laterally  averaged  velocity  inte¬ 
grated  over  h  (m/s) 


Indicator  of  flow  direction  for  upstream 
differencing 

Cell  volume  (B*h‘Ax)(m3) 

Settling  velocity  (m/s) 
z-direction  velocity  (m/s) 

z-direction,  laterally  averaged  velocity  (m/s) 

(W  in  FORTRAN  code) 

Wind  speed  (m/s) 

Cartesian  coordinates:  x  is  along  the  lake  center 
line  at  the  water  surface,  positive  to  the  right, 
and  z  is  positive  downward  from  the  x-axis  (m) 

Horizontal  spatial  step  (m)  (DLX  in  FORTRAN  code) 

Reservoir  bottom  elevation  (m) 

Surface  elevation  (m) 

Density  (kg/m3) 
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